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LONG-TERM  GOALS 

The  goal  is  to  assess  the  ability  of  LES  with  modern  subgrid  turbulence  closures  to  predict  and 
understand  transport  and  mixing  in  stratified  flows  over  boundaries. 


OBJECTIVES 


•  perform  large  eddy  simulation  (LES)  of  canonical  turbulent  flows  over  boundaries  and 
compare  against  available  laboratory,  field  and  direct  numerical  simulation  (DNS)  data 
using  metrics  of  mean  flow  profiles,  Reynolds  stresses  and  dissipation  rate. 

•  quantify  the  effect  of  grid  resolution,  Reynolds  number  and  stratification  level  on  model 
performance. 

•  represent  oceanic  small-scale  mixing  processes  in  the  AESOP  experiment  using 
fine-resolution  LES  and  help  understand  field  data. 


APPROACH 

A  non-hydrostatic  code  that  numerically  solves  the  unsteady,  three-dimensional,  primitive 
equations  is  used.  Advanced  models  such  as  the  dynamic  mixed  model  and  the  dynamic  eddy 
viscosity  model  are  utilized  to  represent  subgrid  processes  in  the  LES  approach.  A  novel  near-wall 
model  has  been  developed  so  as  to  increase  the  Reynolds  number  of  boundary  flows  to 
realistically  large  geophysical  values.  The  code  is  capable  of  capturing  the  effect  of  bottom  slopes, 
small-scale  roughness,  and  oscillatory  forcing. 


WORK  COMPLETED 

Near  wall  model 

When  simulating  geophysical  boundary  layers  at  a  large  Reynolds  number,  in  addition  to  the 
subgrid  model  for  the  outer  layer,  a  near-wall  model  (NWM)  is  inevitable  for  handling  the  inner 
layer.  This  is  due  to  the  fact  that  the  filter  scale  is  typically  much  larger  than  the 
viscous/roughness  scales,  the  generation  of  turbulence  at  the  rough  wall  is  not  resolved,  and 
consequently  the  near-wall  Reynolds  shear  stress  must  be  explicitly  accounted  for  through  a 
model.  This  issue  has  ben  faced  by  the  atmospheric  boundary  layer  (ABL)  community,  where  it  is 
known  that  LES  overpredicts  the  mean  shear  near  the  surface  [7,  5]  relative  to  expected  similarity 
theory.  Procedures  to  circumvent  the  problem  with  a  near- wall  mocenteringdel  have  been 
advanced  [7,  11,  10],  but  there  is  no  well-accepted  solution. 
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Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 
VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 
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Figure  1:  (a)  Mean  velocity  profiles  a 
velocity  for  smooth  wall  channel  flow 
near-wall  model  (NWM)  is  asse 


%d  (b)  streamwise  and  wall-normal  r.m.s 
with  Re*  =  2000.  The  performance  of  the 
ised  by  comparison  with  DNS  data. 


Near-wall  models  for  geophysical  boundary  layers  are  commonly  formulated  by  applying  an 
approximate  boundary  condition  at  the  wall,  an  unsteady  generalization  of  the  log  law.  That 
proposed  by  Marusic  et  al.  [6]  which  will  be  referred  to  here  as  the  MKP  model  has  been 
adopted.  The  MKP  model  is  formulated  as  the  following  approximate  wall  boundary  condition: 

Tw,i(x,y,t)  _  Uj(  1) _ [uj(x  +  As,y,z(l))  -  Uj(  1)] 

p0ul  (17(1)2  +  F(l)2)l/2  P  U*  ’  U 

where  z(  1)  is  the  location  of  the  gridpoint  nearest  the  wall  and  U(l)  and  V(l)  are  the  plane 
averaged  horizontal  velocity  components  at  the  first  gridpoint.  Marusic  et  al.  [6]  hypothesized 
that  the  constant  /3  should  be  universal,  and  empirically  determined  /3  —  0.10.  The  instantaneous 
friction  velocity  and  mean  wall  stress  are  estimated  by  assuming  that  the  plane-averaged  velocity 
at  the  first  gridpoint  follows  the  expected  log  law.  Note  that  the  friction  velocity,  tt*,  and  surface 
stress  angle,  a,  do  not  have  to  be  prescribed  a  priori  but  are  a  product  of  the  simulation. 


The  MKP  near-wall  model  along  with  the  dynamic  eddy  viscosity  model  in  the  exterior  has  been 
tested  by  us  for  smooth  and  rough  channel  flow,  and  for  a  bottom  Ekrnan  layer  with  and  without 
stratification.  In  all  cases,  the  mean  velocity  is  overpredicted  with  respect  to  the  log  law.  This 
trend  is  illustrated  with  channel  flow  at  the  friction  Reynolds  number,  Re*  =  2000,  where  the 
direct  numerical  simulation  (DNS)  of  Hoyas  and  Jimenez  [4]  is  available  for  comparison.  Figure 
(a)  shows  the  mean  velocity  profiles  as  a  function  of  z+  =  zu*/v  for  smooth  wall  channel  flow 
simulations.  The  mean  velocity  from  the  DNS  follows  the  expected  log  law  between  z+  >  30  and 
z/5  <  0.2.  It  is  apparent  from  Figure  (a)  that  the  NWM-LES  without  stochastic  forcing  (blue  line 
with  filled  circles)  overpredicts  the  mean  shear  where  a  log  law  is  expected.  A  novel  solution 
based  on  adaptive  stochastic  forcing  has  been  developed  as  described  in  detail  and  validated  by 
Taylor  and  Sarkar  [15].  A  summary  is  given  below. 


Adaptive  stochastic  forcing 

A  stochastic  forcing  term  is  added  to  the  r.h.s  of  the  wall  normal  momentum  equation  so  as  to 
enhance  the  vertical  velocity  fluctuations  while  modifying  the  local  Reynolds  stress  according  to 
the  expected  value  from  the  log  law.  The  characteristic  length  and  time  scales  of  the  forcing  are 
set  by  the  grid-size  and  time-step,  respectively.  The  proposed  forcing  function  can  be  written 

fz(x,y,z)  =  ±UA(z),  (2) 
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Figure  2:  (a)  Mean  velocity  profiles  for  rough  wall  channel  flow  with  zo/5  =  2.8 1.x  10“ 
and  Re*  =  5600,  (b)  Mean  velocity  profiles  for  a  bottom  Ekman  layer  with 

zq/5  =  5.7  *  10~5  and  Re*  =  60,  00. 


where  TZ  is  a  random  number  between  0  and  1  and  A(z)  is  an  amplitude  function.  The  forcing 
amplitude  is  evaluated  dynamically  in  order  to  ensure  that  the  mean  velocity  follows  the 
logarithmic  law  near  the  wall.  The  difference  between  the  resolved  shear  and  the  expected 
logarithmic  law  value  is 

d(u)2  d{v)2 

dz  dz 


(3) 


The  forcing  amplitude  can  then  be  locally  adjusted  using  a  proportional  controller  with 


A(z)n+1  =  A(z)n  + 


(4) 


where  t  sets  the  relaxation  time  and  should  be  large  enough  to  allow  the  flow  to  adjust  to  the 
forcing,  and  typically  we  have  used  r  ~  5/u*.  The  sign  of  the  r.h.s  of  Eq.  (2)  is  chosen  to  increase 
or  decrease  the  u',  w’  correlation  according  to  the  sign  of  e(z).  The  stochastic  forcing  must  satisfy 
the  following  constraints:  its  amplitude  must  be  small  and  its  width  in  physical  space  be  thin  so 
that  the  “large  eddies”  external  to  the  near-wall  region  are  not  artificially  forced.  These 
constraints  are  met,  for  example,  in  the  rough  Ekman  layer,  the  amplitude  is  less  than  6%  of  the 
r.rn.s  vertical  velocity  and  is  non-zero  only  for  5  points  near  the  boundary,  less  than  8%  of  the 
boundary  layer  thickness. 


Returning  to  the  example  in  Fig.  (a),  it  can  be  seen  that  a  much  better  agreement  with  the  DNS 
in  the  logarithmic  region  is  achieved  when  the  stochastic  forcing  term  is  included.  The  streamwise 
and  wall- normal  turbulent  velocities  are  shown  in  Figure  (b).  When  the  NWM-LES  is  used,  the 
magnitude  of  the  wall-normal  rms  velocity  is  significantly  less  than  the  DNS.  Adding  stochastic 
forcing  increases  the  wall-normal  velocity  fluctuations  as  expected,  and  results  in  a  better 
agreement  with  the  DNS.  The  same  conclusion  is  true  for  the  spanwise  velocity  fluctuations  (not 
shown  for  clarity).  Without  forcing,  the  streamwise  velocity  fluctuations  compare  reasonably  well 
in  magnitude  to  the  DNS,  but  an  improvement  is  still  seen  when  forcing  is  added,  resulting  in  a 
shift  in  the  location  of  the  maximum  near  the  wall  and  an  increase  in  the  outer  layer  fluctuations. 
Fig.  (a)  shows  the  example  of  a  rough  boundary  layer  where  a  comparison  is  made  with 
experiments  of  Bakken  et  al.  [2]  where  wire  mesh  was  used  to  roughen  the  ceiling  and  floor  of  a 
wind  tunnel.  The  friction  Reynolds  number  is  5600,  and  the  roughness  length  scale  is  k+  =  187 
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Table  1:  Physical  Parameters 


Ri* 

Noo/f 

Re* 

Pr 

/Uqq 

zo/5 

«o 

0 

0 

1.08  *  106 

0.0488 

4.80  *  10~5 

15.4° 

1000 

31.6 

1.09  *  106 

5-10 

0.0490 

4.78  *  10~5 

18.9° 

5625 

75 

1.12  *  106 

0.0497 

4.71  *  10-5 

24.8° 

indicating  that  the  flow  is  fully  rough.  Fig.  (b)  shows  an  application  to  an  unstratified  bottom 
Ekrnan  layer  where  a  rough- wall  log- law  is  expected  to  hold  throughout.  In  both  examples,  the 
results  with  forcing  (green  line  with  unfilled  circles)  are  much  better  than  without  (blue  line  with 
filled  circles). 


Bottom  Ekman  layer  under  a  current 

A  stratified  bottom  Ekman  layer  over  a  non-sloping,  rough  surface  has  been  studied  using 
large-eddy  simulation  in  order  to  examine  the  effects  of  the  outer  layer  stratification  on  the 
boundary  layer  structure.  Three-dimensional,  unsteady  simulations  of  a  stratified  bottom  Ekman 
layer  have  not  been  performed  before  although  the  ABL  analog  of  a  boundary  layer  stratified 
through  a  surface  cooling  heat  flux  has  been.  Details  of  our  LES  investigation  are  given  by  Taylor 
and  Sarkar  [14] ;  a  summary  follows  below. 

The  turbulent  Ekman  layer  considered  here  is  formed  when  a  steady  flow  in  geostrophic  balance 
encounters  a  flat  boundary.  The  x-axis  is  aligned  with  the  external,  steady  zonal  flow  in 
geostrophic  balance.  The  seafloor  is  represented  by  a  non-sloping,  rough  surface.  The  roughness 
elements  are  too  small  to  resolve  directly  by  the  grid,  but  their  effect  is  parameterized  through  a 
near-wall  model.  The  seafloor  is  also  assumed  to  be  adiabatic.  The  boundary  conditions  are 
periodic  in  the  horizontal  directions;  a  sponge  layer  along  with  a  non-reflecting  boundary 
condition  is  used  in  the  vertical.  The  domain  size  is  approximately  25  x  25  x  25  where 
5  =  u*/ /  «  100m.  The  parameters  considered  in  this  study  are  listed  in  Table  1.  Note  that  Re oo 
is  the  same  for  each  simulation,  but  the  drag  coefficient  depends  on  the  outer  layer  stratification, 
and  hence  the  friction  Reynolds  number  varies  between  each  case.  We  have  performed  simulations 
at  three  different  values  of  Ri* ,  equivalent  to  changing  the  free- stream  density  gradient.  For 
comparison  with  oceanographic  conditions,  observations  of  the  bottom  boundary  layer  over  the 
Oregon  shelf  by  Perlin  et  al.  [8,  9]  provide  estimates  of  Re*  =  4  *  104  —  8  *  105  and  N/  f  =  75. 
Therefore,  both  the  Reynolds  numbers  and  stratification  levels  considered  in  the  present  study  are 
comparable  with  field  data.  In  order  to  compare  our  simulation  results  to  realistic  conditions,  we 
will  use  Uoo  =  0.2 m/s,  f  =  10-4s“l,  and  v  =  10“6m2/s  which,  assuming  that  u^/Uoo  =  0.05, 
yields  5  ~  100m  and  zq  =  0.5cm.  The  roughness  lengthscale  is  consistent  with  observations  by 
Perlin  et  al.  [8]  who  found  that  zo  =  0.05cm  -  2  cm  depending  on  the  method  used  to  infer  the 
wall  stress. 

The  subgrid  scale  stress  tensor,  r,  and  the  subgrid  scale  density  flux,  A,  are  evaluated  using  the 
dynamic  Smagorinsky  model  [3] .  The  dynamic  procedure  avoids  the  empirical  specification  of  the 
Smagoinsky  coefficient  and  has  been  shown  to  perform  well  for  wall-bounded  flows  and  density 
stratified  flows  [1,  12].  Details  of  the  basic  LES  model  can  be  found  in  Armenio  and  Sarkar  [1]. 

To  avoid  the  need  to  resolve  the  very  small  turbulent  motions  near  the  lower  wall,  we  have  used 
the  near- wall  model  (NWM)  described  previously  in  this  report. 

Mean  properties  of  the  boundary  layer 

The  plane-averaged  temperature  profiles  after  the  flow  has  reached  quasi-steady  state  are  shown 
in  Figure  3(a).  The  components  of  the  Reynolds  averaged  horizontal  velocity  are  shown  in  Figure 
3(b).  Comparing  with  the  temperature  profiles,  it  is  apparent  that  most  of  the  Ekman  transport 
is  confined  to  the  mixed  layer.  The  increase  in  cross  stream  velocity  results  in  a  broadening  of  the 
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Figure  3:  (a)  Plane  averaged  temperature  profiles,  (b)  Plane  and  time-averaged 
horizontal  velocity  components,  (c)  Angle  of  Ekman  veering.  Note  that  S  &  100  m  in 

this  and  the  following  figures. 


Ekman  spiral  in  a  hogograph  (not  shown).  It  is  interesting  to  note  that  while  the  density  gradient 
is  zero  near  the  lower  wall,  the  surface  turning  angle  ao  =  tan~l(Ty/rx )  increases  with  the  outer 
layer  stratification,  as  listed  in  Table  1.  The  angle  of  Ekman  veering,  defined  by  a  =  tan -1 
is  shown  as  a  function  of  z/5  in  Figure  3(c).  When  the  flow  is  unstratified,  the  turning  of  the 
mean  velocity  occurs  gradually  throughout  the  boundary  layer.  The  rate  of  turning  in  the  mixed 
layer  does  not  depend  strongly  on  the  outer  layer  stratification.  Therefore,  since  the  magnitude  of 
the  turning  angle  at  the  wall  increases  with  the  outer  layer  stratification,  and  the  boundary  layer 
height  decreases  significantly,  the  rate  of  turning,  da/dz,  becomes  very  large  in  the  pycnocline. 


The  mixed  layer  height  decreases  strongly  with  increasing  N  with  a  value  of  15  m  at  N/ f  =  75, 
Fig.  3(a).  The  boundary  layer  height  based  on  the  mean  velocity  also  decreases  strongly  with 
increasing  stratification  as  can  be  inferred  from  Fig.  3(b).  Zilitinevich  and  co-workers  have 
proposed  a  scaling  law  based  on  a  combination  of  theory,  field  data  [16]  and  simulations  [17]: 


hzE 


„  Ujf  (  CftCuN  -/Vqc 

Cr7  1 1  +  lF7 


(5) 


where  the  constants  Cr  =  0.5  and  Cun/C |  =  0.56  were  determined  by  fitting  to  data  from 
large-eddy  simulation.  We  define  the  boundary  layer  height  as  the  location  where  the  Reynolds 

stress  ( (u'w1)2  +  (v'w1)2)  =  0.1  uf  This  estimate  agrees  well  with  the  values  given  by  the 


correlation,  Eq.  (5),  namely,  a  boundary  layer  depth  of  h/5  =  0.215  and  h/5  =  0.147  for 
N/f  =  31.6  and  N/f  =  75. 


The  density  gradient  intensifies  at  the  top  of  the  mixed  layer  to  form  a  bottom  pycnocline,  shown 
in  Figure  4(a).  The  zonal  velocity  has  the  following  interesting  feature:  a  zonal  jet  with  higher 
than  the  outer  geostrophic  velocity  forms  in  the  pycnocline,  Figure  4(b).  A  peculiar  feature  of  the 
mean  velocity  profile  when  N/f  =  75  is  that  the  mean  velocity  and  the  mean  shear  are  maximum 
at  the  same  location  near  the  center  of  the  pycnocline,  Figure  4(c). 


5 


(a) 


(b) 


(0 


- A///'=31.6  - N/f=15 


Figure  4:  (a)  Temperature  gradient,  (b)  horizontal  velocity  magnitude  (c)  mean  shear 


(a) 


Figure  5:  Reynolds  averaged  horizontal  velocity  magnitude 
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The  law-of-the-wall  is  expected  to  hold  in  a  region  far  enough  from  the  wall  where 
viscosity/roughness  scale  can  be  neglected,  but  near  enough  to  the  wall  so  that  the  boundary  layer 
depth  is  not  felt.  The  Reynolds  averaged  horizontal  velocity  magnitude  on  a  semi-logarithmic 
scale  is  shown  in  Figure  5(a).  Very  near  the  wall,  all  cases  are  in  reasonably  good  agreement  with 
the  unstratified  logarithmic  law.  Deviations  from  the  logarithmic  velocity  profile  in  the  cases 
when  stratification  is  present  can  be  seen  clearly  by  calculating  the  normalized  velocity  gradient, 


$  = 


KZ 

U* 


d  <  u  > 
dz 


(6) 


The  quantity  $  can  be  interpreted  as  the  ratio  of  the  observed  mean  velocity  gradient  to  that 
expected  from  the  logarithmic  law.  When  the  outer  layer  is  stratified,  the  mean  shear  in  the 
pycnocline  increases  significantly  to  $  =  2.  It  is  worth  noting  that  deviations  from  the 
law-of-the-wall  begin  well  within  the  mixed  layer. 


Turbulence  properties  of  the  boundary  layer 


One  consequence  of  the  damping  of  turbulence  by  stratification  is  a  decrease  in  the  turbulent 
stresses,  <  v'w'  >  and  <  v'w'  >.  The  decrease  in  boundary  layer  height  with  increasing 
stratification  is  very  apparent  from  the  Reynolds  stress  profiles  shown  in  Figure  5(b).  Here,  the 
Reynolds  stress  includes  both  the  resolved  and  subgrid  scale  contributions.  Note  that  above  the 
boundary  layer,  the  Reynolds  stress  approaches  a  small,  nonzero  value  owing  to  the  presence  of 
turbulence-generated  internal  gravity  waves.  However,  the  Reynolds  stress  associated  with  the 
waves  is  much  smaller  than  in  the  turbulent  region,  so  this  does  not  interfere  with  the  definition 
of  boundary  layer  height.  It  is  evident  from  Figure  5(b)  that  the  Reynolds  stresses  change  more 
rapidly  with  height  when  the  outer  layer  stratification  is  stronger.  It  can  be  shown  that  the 
behavior  of  the  Reynolds  stresses  in  the  outer  layer  is  consistent  with  the  formation  of  the  zonal 
jet  observed  in  Fig.  4(b). 

The  plane-averaged  turbulent  kinetic  energy  budget  can  be  found  by  dotting  lT  into  the 
perturbation  momentum  equations: 


1  8  <  u'ju'i  > 

2  dt 


Id  ,  ,  ,  d  ,  , 

“2 dz<WUiUi>-dz<WV  >  + 


i  d2k  ^  a 
r^Yz 2“  <  ij 


.  .  /  /  . 

><  UtUj  >  — 


1 


TT  /  /  d  ,  lyu,.'  i jru 

-Rlr  <W  p  >  <  UiT3 1  >  +  <  Tji—  >=  —  =  0. 


Re 
du ' 


du'  du' 

<  — - — -  > 
dxj  dxj 


dxj 


dk 

dt 


(7) 


Reading  from  left  to  right,  the  terms  on  the  right  hand  side  of  Eq.  (7)  can  be  identified  as  the 
turbulent  transport,  pressure  transport,  viscous  diffusion,  production,  dissipation,  buoyancy  flux, 
subgrid  transport  and  subgrid  dissipation,  respectively.  The  leading  terms  in  the  turbulent  kinetic 
energy  budget  for  N^/f  =  75  is  shown  in  Figure  .  Near  the  wall,  the  leading  order  balance  is 
between  production  and  dissipation  and  does  not  differ  significantly  from  the  unstratified  case  as 
shown  in  the  inset.  In  the  upper  portion  of  the  mixed  layer,  the  turbulent  transport  appears  as  a 
source  term  representing  the  advection  of  turbulent  eddies  towards  the  pycnocline.  Pressure 
transport,  buoyancy  flux,  and  dissipation  represent  the  dominant  energy  sinks  in  the  upper  mixed 
layer.  In  the  pycnocline,  starting  at  about  z/5  =  0.15,  the  turbulent  transport  and  buoyancy  flux 
decrease  as  stratification  suppresses  the  turbulent  motion  while  pressure  transport  becomes  the 
dominant  source  term.  When  the  pressure  transport  is  positive  the  vertical  energy  flux,  <  p'w'  > 
is  increasing  with  height,  consistent  with  an  internal  wave  field  that  is  gaining  energy.  This  is 
direct  evidence  of  the  generation  of  internal  waves  by  the  interaction  between  boundary  layer 
turbulence  and  a  stable  stratification,  as  examined  in  detail  by  us  in  a  previous  study  [13]. 
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IMPACT/ APPLICATIONS 
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Figure  6:  Turbulent  kinetic  energy  budget  at  the  top  of  the  mixed  layer  and 
pycnocline  for  Noo/f  =  75.  Inset  shows  the  TKE  budget  near  the  wall. 


The  detailed  mean  and  turbulence  structure  of  a  stratified  BBL  has  been  obtained  using  LES. 
The  height  of  the  BBL  and  the  drag  coefficient  are  outputs  and,  thus,  can  be  used  as  inputs  for 
regional  models.  A  new  near-wall  model  with  adaptive  stochastic  forcing  has  been  developed  to 
solve  the  problem  of  mean  velocity  over-prediction  near  the  boundary  in  previous  simulations. 
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